APS/123-QED 



Spontaneous Magnetic Ordering in a Ferromagnetic Spinor Dipolar Bose-Einstein 
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We study the spin dynamics in a spin-1 ferromagnetic Bose-Einstein condensate with magnetic 
dipole-dipole interaction (MDDI) based on the Gross-Pitaevskii and Bogoliubov theories. We find 
that various magnetic structures such as checkerboards and stripes emerge in the course of the 
dynamics due to the combined effects of spin-exchange interaction, MDDI, quadratic Zeeman and 
finite-size effects, and non-stationary initial conditions. However, the short-range magnetic order 
observed by the Berkeley group [Phys. Rev. Lett. 100, 170403 (2008)] is not fully reproduced 
in our calculations; the periodicity of the order differs by a factor of three and the checkerboard 
pattern eventually dissolves in our numerical simulations. Possible reasons for the discrepancy are 
discussed. 

PACS numbers: 03.75.Mn,03.75.Kk,67.30.he 



I. INTRODUCTION 

Experiments on dilute Bose-Einstein condensates 
(BECs) have exhibited a rich variety of phenomena, most 
of which have successfully been explained by theory, in- 
cluding the system with a strong magnetic dipole-dipole 
interaction (MDDI) [J. The magnetic crystallization of a 
spin-1 ferromagnetic 87 Rb condensate recently observed 
by the Berkeley group @, Q provides one of the few 
anomalies that have so far defied theoretical explanation. 
It has been argued that the MDDI, which is long-ranged 
and anisotropic, plays a pivotal role in the magnetic crys- 
tallization [2|-|8j]; however, no satisfactory account of the 
experiment has been presented. 

The aim of this paper is to clarify what the mean- 
field and Bogoliubov theories predict under the condi- 
tions of the Berkeley experiment 0- The authors of 
Ref. have found that a helical spin structure of a quasi- 
two-dimensional (2D) condensate spontaneously develops 
into a short-range crystalline pattern of magnetic do- 
mains in a time scale of a few hundreds of milliseconds. 
The typical size of the magnetic pattern is A CX p — 10 /jm 
(the size of the magnetic domain is A CX p/2), and the 
growth rate of the crystalline pattern increases with de- 
creasing the pitch of the initial spin helix. On the other 
hand, when the initial spin configuration is uniform, the 
short-range magnetic patterns do not emerge over a pe- 
riod of 300 ms. During the dynamics, the condensate 
remains fully magnetized and the longitudinal magne- 
tization is much smaller than the transverse magneti- 
zation. They also found that the number of spin vor- 
tices increases as the short-range magnetic patterns de- 
velop. The nature of the pattern is insensitive to the 
strength of the quadratic Zeeman energy q in the range 
of 0.8 < q/h < 4 Hz. The same group also reported 



the emergence of a similar crystalline pattern by cool- 
ing initially unmagnetized gases [||. This result suggests 
that the system of spin-1 87 Rb BEC possesses an intrinsic 
mechanism that stabilizes the crystalline magnetic order 
with a characteristic size of A exp . 

The complexity of this system arises from the inter- 
play among three effects related to the magnetism: the 
short-range ferromagnetic spin-exchange interaction, the 
quadratic Zeeman effect which favors transverse magne- 
tization [ij [l(|, and the long-range anisotropic MDDI 
which induces spatial spin textures [ll]-[l3|. In addition 
to them, the initial conditions and noises, the finite-size 
effect and nonuniform density profile due to the trapping 
potential contribute to the dynamics. In this paper, we 
take into account all these features of the spinor dipolar 
BEC and investigate the extent to which we can under- 
stand the observed phenomena. We do find the emer- 
gence of magnetic checkerboard patterns due to the in- 
homogeneity and finite-size effects of the trapping poten- 
tial; however, the periodicity differs by at least a factor of 
three and the checkerboard pattern eventually dissolves 
in our numerical simulations. We find that the MDDI in- 
duces magnetic patterns that are not checkerboard pat- 
terns but spin-helix or staggered-domain structures de- 
pending on the magnitude of the quadratic Zeeman en- 
ergy 

This paper is organized as follows. In Sec. [Til we de- 
scribe a system of the spin-1 spinor dipolar BEC, and 
review the ground-state properties in the absence of the 
MDDI. In Sec. IIII1 we investigate the linear stability of 
the spin-1 spinor dipolar BEC in an infinite quasi-2D 
system. We analytically solve the Bogoliubov equation 
for a uniform spin structure. The Bogoliubov spectrum 
for a helical spin structure is numerically obtained. In 
Sec. IIV1 we discuss the spin dynamics in an oblate trap. 
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We numerically solve the three-dimensional (3D) Gross- 
Pitaevskii equation (GPE), and examine the individual 
effect of the trapping potential, initial spin helix, and 
the MDDI. In Sec. [V] we compare the obtained results 
with the Berkeley experiment (2|. The possible reasons 
for the discrepancy between the experimental results and 
our results are discussed in Sec. IVI1 We provide conclud- 
ing remarks in Sec. IVIII 



where ui^ = gp^B/h with gp being the hypcrfme g- 
factor and /xb the Bohr magneton. In this paper, we take 
the spin quantization axis z along the external magnetic 
field, i.e., es = z. Then, the linear Zeeman energy is 
expressed as 

Ezi = /"dr^fiw L m|* m (r)| 2 . (8) 



II. SPIN-1 SPINOR DIPOLAR BEC 
A. Mean-field energy 

We consider a spin-1 BEC of N atoms confined in an 
optical trap t/trap(f")- The zero-temperature mean-field 
energy is given by 



E= fdr ]T *^(r) 

>> — — n _1_ i 



m=0,±l 

E s + Ezi + Ez2 



h 2 v 2 

2M 

Edd 



+ ^trap(r) 



* m (r) 



(1) 



where M is the atomic mass, \l/ m (r) is the order parame- 
ter of the condensate with magnetic sublevel m = 0,±1, 
and E s , Ezi, Ez2 and -Edd are the short-range interaction 
energy, the linear and quadratic Zccman energy, and the 
long-range MDDI energy, respectively. The order param- 
eter is normalized to satisfy J2 m I dr\^ m \ 2 = N. 
The short-range interaction energy is given by 



E - = \J dr [co^ 2 (r-)+ Cl |/(r)| 2 ] , 



where 



n(r) = l*m(r) 

nt 

is the atom-number density and 



f(r) 



V* m (r)F mm >V m >(r) 



(2) 



(3) 



(4) 



is the spin density with F = (F x , F y ,F z ) being the vector 
of the spin-1 matrices, and the interaction coefficients in 
Eq. © are given by 



Co 

Cl 



Anh 2 ao + 2o2 
~~M 3 : 
Airti 2 ci2 — ciq 
~M 3 ' 



(5) 
(6) 



with as (S = 0, 2) being the s-wave scattering length for 
the scattering channel with total spin S. 

In the presence of an external magnetic field B = Be b , 
the linear Zeeman energy is given by 

Ezi = f dr Yl ' F )™»' *™'( r )< ( ? ) 



The quadratic Zccman energy is induced by a linearly 
polarized microwave field as well as by an external mag- 
netic field as H H 



Ez2= J drJ2** m (r) 



q B (e B ■ F) 2 + q EM (e EM • F) 2 * m /(r), (9) 



where i-EM is the direction of the polarization of the mi- 
crowave field, qs = (sfMb-B) 2 / -E^f with E^f being the 
hyperfine splitting energy, and <7em = — h 2 fl 2 /(4S) with 
f2 being the Rabi frequency and 8 the detuning. In this 
paper, we take eEM to be parallel to es = z. The 
quadratic Zeeman energy then becomes 

Ez2= f drJ2im 2 \^ m (r)\ 2 , (10) 

m 

where q = q B + qEM- 

The general form of the MDDI energy is given by 

E dd =c dd fdr I ' dr'Q vv ,{r-r')f v {r)f v ,{r'), (11) 



where c dd = Mo(5fA*b) 2 /(47r) with /io being the magnetic 
permeability of the vacuum, and Q UV '{r) is the dipole 
kernel which will be given in Sec. Ill CI Here and hence- 
forth, the Greek subscripts that appear twice are to be 
summed over x,y, and z, 

Note that if the atomic cloud is isolated in the vac- 
uum, the total magnetization along the external mag- 
netic field is conserved, as long as the dipolar relaxation 
(the spin relaxation due to the MDDI) can be ignored. 
The dipolar relaxation is dominant for atoms with large 
magnetic dipole moments such as 52 Cr atoms [lj|, while 
it is negligible in BECs of alkali atoms [ItJ. In the lat- 
ter case, the linear Zeeman term can be eliminated if 
we choose the rotating frame of reference in spin space 
with the Larmor frequency wl by transforming \& m (r, t) 
to * m (r,i)e" lmtJL *. In the absence of the MDDI, the 
total energy functional (TTJ) is invariant under this trans- 
formation since E s and Ez2 are invariant. However, in 
the presence of the MDDI, Edd is not invariant under this 
transformation, and therefore we have to use the mod- 
ified dipole kernel in the rotating frame, which will be 
given in Sec. Ill CI 
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B. Nonlocal Gross-Pitaevskii equation 

The mean-field dynamics of the system is governed by 
the nonlocal GPE: 

..d T S(E — fj,N) 
dt m SV* 



2M 



V 2 + C/ trap (r) - n + qm 2 + c n 



■nun tfm', (12) 



m'=0,±l 

where /i is the chemical potential and 

K{r) = J dr'Q^ir-r 1 ).^^ 1 ) 

is the dipole field. 

We consider a harmonic trap 



(13) 



Utrap(r) = M {u\x\ + u 2 2 x\ + oj^xl) /2, (14) 

where xt = ■ r with e.% being the trap axis. Assum- 
ing u 3 3> Wi,2j the wave function in the e 3 direction is 
approximated by a Gaussian 



h(x 3 ) 



1 



(2^rf 2 )i/4 Cxp y 4d 2 , 

and the order parameter can be written as 

$m(r,t) =tp m (x 1 ,X2,t)h(x 3 ). (16) 

Multiplying Eq. (fl"2"|) by h(x 3 ) and integrating over x 3 , 
we obtain the quasi-2D GPE: 

at 



(15) 



(17) 



m'=0,±l 



^J) ^trap' > ( T *-L) 



where rj_ = (xi,^), V 2 
4f (iofxf + uj 2 x 2 ) , c ,i = c ,i/V47rd 2 and Cdd = 

Cdd/V 47rc? 2 . We define the 2D number density, spin den- 
sity, and dipole field by 



dx 3 n(r), 
dx 3 f{r), 



2„' n( 2D 



(18) 
(19) 
(20) 



respectively, where 
Qlv (r-r) 



VAird 2 / / dx 3 dx' 3 h 2 {x 3 )h 2 {x' 3 )Q^(r - r'). (21) 



Using the Fourier transform of the dipole kernels 

Qw<{r)=Y,Qkvwe ikr , (22) 

€ ) (--i)=EC/ ifei - ri . (23) 

fe_L 

with fej^ = (fci, fc 2 ) = ((ei • fc), (e 2 •&)), Eq. (f2Tj) is rewrit- 
ten as 

ft = ^/ * 3 e- d2fc 3Q w , (24) 
where ^3 = <>3 • fe. 

C. Dipole kernel 

Since the MDDI is long-ranged and couples spin and 
orbital degrees of freedom, its kernel depends on the ge- 
ometry of the condensate and on the frame of reference 
in spin space. In this subsection, we discuss the dipole 
kernel in 3D and quasi-2D condensates both in the lab- 
oratory frame and the rotating frame at the Larmor fre- 
quency. 

The dipole kernel in the laboratory frame of reference 
is given by 



"ir v r v > 



(25) 



where r = \r\, f = r/r, and its Fourier transform is given 
by 



A(lab) 
^kvv 1 



(26) 



where fe = fc/|fc|. To calculate the 2D dipole kernel in 
the laboratory frame, we expand fc„ as fc„ = Y^=i ki{e.i) v 
and substitute Eq. (J2l>| in Eq. flUJ). Using the following 
integrals: 



dh 



3e 



-dH--, 



fir J k 2 + k% + k 2 

f G(k±d)(k x )i(k±)j = 1,2), 
= I l-G(fcxd) (i = i = 3), 

[ (otherwise), 

with k± = \k±\, k± = k±/k±, (fej_)j = k± ■ e-i, and 



(27) 



G{k) = 2ke k 



dt, 



(28) 



we obtain 

A(2D,lab) 
^ k± vv' 



4tt 



[5 vv i - 3{e 3 ) v (e.3,) v i 



+ 4irG{k ± d) {k±)„(k ± ) v > - (e 3 )„(e 3 ),. 



(29) 
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where (k±) v = YU=i 2v i -^)i{ l ^i)"- ^ can be sn °wn that 
G(k) is a monotonically increasing function that satisfies 
G(0) = and G(oo) = 1. 

Next, we consider the dipole kernel in the rotating 
frame of reference. Since the spin density vector in the 
rotating frame /( rot ) is related to that in the laboratory 
frame /( lab ) as 



where 



R 



/■(lab) _ d f(rot) 
Jv lx vv' J v i 1 



'coswLt — sinwLi 0' 
sinwL^ coso>l£ 
1, 



(30) 



(31) 



the dipole kernel in the rotating frame of reference is 
given by R„ p Q pp* (r)R p > v > , where T denotes the matrix 
transpose. When the Larmor precession is much faster 
than the dynamics caused by the MDDI, we can use the 
time-avera ged dipole kernel over the period of the Larmor 
precession [la |: 



R T n (lab) 
KpQpp' 



(r)R 



p' v' 



1 1 

2 



(<W - 35 ZU 5 ZU >) , 



(32) 



where (■ • • ) represents the time average over 2tt /ui^. This 
is the case with the Berkeley experiment Q, where the 
MDDI energy is c^n/h ~ 1 Hz and the Larmor fre- 
quency is il1l/(27t) = 115 kHz under the magnetic field 
of B = 165 mG. This approximation ignores the spin- 
orbit coupling terms that induce spin relaxation via the 
Einstein-de Haas effect [l9l - |2~i1 ]. and therefore the spin 
and orbital parts are decoupled in Eq. (|32]l. The ex- 
perimental finding that the longitudinal magnetization 
is conserved over a period of 1 s in 87 Rb BECs [ij} is 
consistent with this approximation. 



The Fourier transform of QJ°?\t) is given by 



j(rot) 



f (l - ikl) (5, 



(33) 



Substituting Eq. (|3"3"|) in Eq. (|2~4")) . we obtain the time- 
averaged 2D dipole kernel in the rotating frame as 



Qk^T ] = (<W - Mz*S av , ) Q k± , (34) 



where 



Q k± = y {l - 3(e 3 )l 3G(k ± d) [(k ± )l - (e 3 )l\ } 



D. 



(35) 



Ground-state phase diagram in the absence of 
MDDI 



We briefly review ground-state properties of the spin-1 
BEC in a uniform system (t/trap = 0) in the absence of 
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FIG. 1: (Color online) Phase diagram of a spin-1 BEC in the 
absence of the MDDI, where the total longitudinal magneti- 
zation is fixed to be zero. 



the MDDI. The phase diagram is shown in Fig. [TJ The 
ground-state spin configuration for the case of q = is 
uniform due to the kinetic term in Eq. ((T|); for c\ > 0, the 
condensate is unmagnetized (/ = 0), while for c\ < 0, 
the condensate is fully magnetized (|/| = n); the former 
phase is called polar or antiferromagnetic, while the latter 
is called ferromagnetic [221 . [23| . 

Though the direction of the spontaneous magnetiza- 
tion in the ferromagnetic phase is arbitrary for q = 0, the 
quadratic Zeeman effect restricts the direction of magne- 
tization 0, [l(| ■ Here we consider the case in which the 
total longitudinal magnetization is fixed to be zero. Sub- 
stituting |*i | = = y/(n - |* | 2 )/2 and minimiz- 
ing Eq. (flj, one finds that the phase transition occurs at 
q = 2|ci|n. For q > 2|ci|n, the quadratic Zeeman energy 
dominates the system and all atoms are condensed in the 
m = state, while the ground state for < q < 2|ci|n 
is partially magnetized in the direction perpendicular to 
the external field. The order parameter for the latter 
case is given by 



(t 




1 



2| Cl |r 



2|ci|n 



V 



1 



(36) 



2|ci|n / 



where </> and a are arbitrary real numbers. The magne- 
tization for this state is given by 



f+ = f w + if y = ne ia \l- 



2\d\n 



(37) 
(38) 



Hence, a denotes the direction of the magnetization in 
the x—y plane, and the amplitude of the magnetiza- 
tion depends on q. Since the spontaneous magnetization 
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breaks the axisymmetry around the external field, this 
phase is called the broken-axisymmetry (BA) phase ■ 
The dynamics of this quantum phase transition has been 
investigated in Refs. [R |2~4T - |28| . When q < 0, the spin- 
polarized state with m — 1 or — 1 can minimize both the 
ferromagnetic interaction and the quadratic Zeeman en- 
ergy. To satisfy the conservation of the total longitudinal 
magnetization, the phase separation of two domains with 
f z = l and —1 must occur. 

The order parameter for the polar phase can be char- 
acterized with a unit vector d as 



J2(F ■ d), 



0. 



(39) 



Since a in Eq. (|36|) specifics the direction of the trans- 
verse magnetization, the helical spin structure can be de- 
scribed as 



tf)(r) = \fn 



in-r 


\ 


Ci(r) 





1 


Co(r) 





e %Kr l 


\C-i(r 



VnKC(r), 



(41) 



where i/> — (ipi, tpo, 4>-i) T , C is a three-component 
spinor satisfying J2 m =o,±i ICm| 2 = 1, and n = 

V2TrcPn(xi,X2,X3 = 0) is assumed to be constant. Sub- 
stituting Eq. dHJ into quasi-2D GPE ([17]). we obtain 



For example, the order parameter for d = z is given 
by (0,1, 0) T . In the absence of the quadratic Zeeman 
effect, the direction of d is arbitrary. However, when 
q is nonzero, the quadratic Zeeman effect restricts the 
direction of d: d = z for q > 0, while d _L z for q < 0. 
The order parameter for the latter case is given by 



„i4> ( d x ~t~ idy 

71 - - 

\ d x + idy 



(40) 



where </> is an arbitrary real number. 

In the presence of the MDDI, the regions of the fer- 
romagnetic phase and the BA phase expand since the 
MDDI favors ferromagnetic structure by developing spin 
textures. The phase diagram in a quasi-2D system with 
a uniform spin structure is investigated in Ref . Q , where 
the phase boundary is predicted to depend on the ge- 
ometry of the system and spin textures, as will also be 
discussed below. 

In the following sections, we consider a BEC of spin-1 
87 Rb atoms which is ferromagnetic with the scattering 
lengths of ao = 101. 8ob and a 2 = 100. 4ae where ob is 
the Bohr radius [29j |. The hyperfine g-factor of the spin-1 
87 Rb atom is gp = —1/2. We choose q to be in the region 
of the BA phase, i.e., < q < 2\ci\n, as in the Berkeley 
experiment 0]. 



III. DYNAMICAL INSTABILITY IN A 
UNIFORM QUASI-2D SYSTEM 



4< 



-—(V- lK F z ) 2 -[l + qF* 



C, (42) 



where 

A k = [ dre~ lk - 



^Qk+^F+OF. + iQ fc - K (C f F_C)F+ 



(43) 



with F + = Fl = F x 
relations: 



iF y , and we used the following 



K^F ± K = e ±tK - r F ± , 
K^F Z K = F Z , 



F 2 = F U F V = 



F+F_ +F_F + 



F 



(44) 
(45) 

(46) 



We consider a spin- helix state with f z = (i.e., 
|Ci | = |C-i I) as in the case of the Berkeley experiment 0]. 
Assuming that C is uniform and satisfies |£i| = |C-i|j we 
can rewrite Eq. (|4"2"]l as 



2M 



F 



c n 



(ci+c dd Q K )n (C f ^C)^ 



u—x,y 



C = K, (47) 



The linear stability of this system has previously been 
discussed in Refs. [1, H, both the initial helix config- 
uration and the MDDI induce the dynamical instability. 
Here we investigate the stability of a spin helix with the 
MDDI in a uniform quasi- 2D system (U^J = 0). For 
the case of a uniform spin structure, our result agrees 
with that in Ref. @. In this section, we consider only 
the quasi-2D system and omit the subscript _L. We as- 
sume that the external magnetic field is much stronger 
than the dipole field and use the dipole kernel (fM)) . 



where we used C^zC = and Q K = <2- K - Note that if 
we replace q with q+H 2 \K\ 2 /(2M) and Si with ci+CddQ/t- 
Eq. g7]) takes the same form as the GPE dTTJ) for the 
stationary state [i.e., dtp m /(dt) = 0] in the absence of 
the MDDI and spin helix. Therefore, as in the case of 
Eq. dSH]), the solution for Eq. (g7J) is given by 



Jtt> 



Co = ^ y/2( l + y , 



(48) 
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where 



q + H 2 \K\ 2 /(2M) 
2n(|ci| - c dd Q K ) 



(49) 



The chemical potential and the energy per particle are 
given by 

q H 2 \k\ 2 

M=2+ AM + h (co +ci + c dd Q K ), (50) 



E K = q - 
2 



1 



+ 2 n ( c ° + Cl + c ddQ K )- (51) 



Note that when the magnetic field is applied parallel 
to the 2D plane (e^ = y), the helical spin structure can 
have lower energy than the uniform spin structure. On 
the other hand, when the magnetic field is perpendicular 
to the 2D plane (e>3 = z), the uniform spin structure 
has the lowest energy. With d = 1.0 /im, n = 2.3 x 
10 14 cm~ 3 , and = y, E K takes a minimum of E Ko — 
E Q = -0.067 Hz x h at k ~ 2tt/(138 (j,m)z. 

The Bogoliubov equation for a spinor BEC is derived 
by substituting 

C = Co + ^[ M fc e4(fe ' r ~ efe * /fi) + vle-'V-r-W)] (52) 



into Eq. (|42|) and linearizing the result with respect to Uk 
and Vk. Here it*, and Vk are three component spinors, 
and therefore, we have the 6x6 eigenvalue matrix equa- 
tion: 



Mi 



N, 



-N*_ k -M*_ k 



Vk) £fc V «fe 



(53) 



where Mk and TV*, are 3x3 matrices defined by 



M k = —(k-KF z y-» + qF 2 



c n 



1 + (Cod) 



(tlF^Fv + F„(CoCl)F„ 



2c dd n 

Cdd"- 



Qk(CIf+Co)F- + q k (CIf-Co)f+ 



-Q k+K F„(CoCl)F+ + Qk- K F + (CoCl)F- 



N k = c n(CoCo T ) + cinF,(CoCo T )^ T 



(54) 



2cd d nSfc^(CoCo T )^ 



Cdd" 



Q fc+K F_(CoCoV- + Q fc - K F + (CoCo T )^+ 



(55) 

Here, (CC*) an d (CC T ) are 3x3 matrices whose (m,m') 
components are given by Cm Cm' anc ^ CmCm'j respectively, 
while C f ^C = Emm' (m(Fv)mrn'(m' is a scalar. 

Since the 6x6 matrix in Eq. (|53|) is not Hcrmitian 
(though Mfe and N k are Hermitian), the eigenvalue can 



be complex. When the eigenvalue has a nonzero imagi- 
nary part, the corresponding mode becomes dynamically 
unstable and exponentially grows or decays. When e& 
is a real eigenvalue with wavenumber fc, — e/~ is also an 
eigenvalue of Eq. (|53[) with wavenumber — fc; the sign 
of 6^ is determined so that the corresponding cigenstate 
satisfies u k Uk — v k Vk = 1. In the presence of the energy 
dissipation, the mode for < is energetically unstable 
due to the Landau instability. 

When k = 0, one of the Bogoliubov modes can be 
obtained analytically. By rewriting u and v as 



(56) 
(57) 



the (u x , v x ) mode is decoupled from the other two modes, 
and the eigenvalue equation (|53| is reduced to the 2x2 
matrix equation given by 





fk 

~9k 



9k 

~fk 



^xk 
V x k 



<k 



V'xk 
Vxk 



(58) 



where 



"2M 2^ idT1 2 ' ' 



i - I — . - _-3 + go a 
9k = \ci\nqo + c dd n Qk- 



The eigenvalue of this spin mode is obtained as 
4 =(fk + 9k)(fk - 9k) 



(60) 



h 2 k 2 

2M + 

h 2 k 2 

2M 



q - c dd n(l - q )(2Qk + So 
+ c dd n(l + qo)(Q k - Qo) 



(61) 



which agrees with the result obtained in Ref. (f|. 

For k ^ 0, we numerically solve the Bogoliubov equa- 
tion in Eq. ([53]) . The results for the case of £3 = y and 
k = (2n/\)z are summarized in Fig. [21 where we show 
the distributions of \lmek/h\ (red) and — Reek/h (blue) 
in the momentum space which correspond to the dynami- 
cal instability and Landau instability, respectively. Since 
the eigenvalue equation (|53|) is the 6x6 matrix equa- 
tion, there are three independent solutions. We have 
numerically confirmed that the lowest-energy mode be- 
comes unstable in some parameter regime (see Fig. [2]), 
and the other two modes are stable. The eigenvalue of 
the lowest-energy mode continuously approaches Eq. (|61[) 
asA-> 00. 

The magnetic fluctuations are analyzed as follows. The 
cigenmodcu; = tt fc e^ fc ' T ' _efct / r ^ +vte~ % ( k ' r ~' L k t / h " 1 induces 
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( a) 9=0.46 Hz, A=30^ m (b) g=1.8Hz. A=30fim (c) g=2.9 Hz. A=30 jim 
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(g) g=0.46 Hz, A=120^m (li) g=1.8Hz, A=120 fim (i) 9=2.9 Hz, A=120 /iia 



029 s 1 



(k) g=1.8Hz,A=3C. (1) g=2.9 Hz, A=oc 



!o.44 s' 1 



t 

054 s 1 



-0.05 0.0 0.05 
k x /2TT [|im J ] 

FIG. 2: (Color) Quadratic Zeeman energy q and helical pitch 
A dependence of the dynamical and Landau instabilities. The 
abscissa and ordinate show k x /2n and k z /2n, respectively. 
Shown are |Imefc| (red) and — Ree*. (blue) which correspond 
to the dynamical and Landau instabilities, respectively; they 
are calculated for (ei, e2, £3) = (z, x, y), k = (2ir/A)z, n(x3 = 
0) = 2.3 x 10 14 cm -3 and d = 1.0 lira. No instability is found 
for (c), (f), and (i). 



fluctuations in the transverse and longitudinal magneti- 
zations as 



A4 = (Co + ™) t ^ y (Co+™) 

— -\Jl + q K Im(wi - w-i) 
Af z = (Co + w)^F z (Co + w) 



— \/l - q K Re(wi - w-i), 
which can be rewritten in the form of 

Afy,z = A yjZ sm[k ■ r - (Re e k )t/h + x y ,z\- 



(62) 
(63) 

(64) 



We may use this result to define the direction of the mag- 
netic fluctuation 9 as 



arctan ■ 



A, 



(65) 



We will call the magnetic fluctuation transverse if < 
9 < 7r/4 and longitudinal if 7r/4 < 9 < n/2. For an 




-0.05 0.0 0.05 
fcr/27r [^nr 1 ] 



I 



tt/2 



transverse Q longitudinal 



FIG. 3: (Color) Direction 9 of the magnetic fluctuation de- 
fined in Eq. (|65|) for the lowest-energy modes. Each panel 
corresponds to that in Fig. [2] The parameters are the same 
as those in Fig. [2] 



eigenmode with a real eigenvalue, u and v are also real, 
and 9 is calculated as 



tan£> 



1 - q K 



u\ + v\ — U-i — V-\ 



u\ — vi — «_i + V-i 



(66) 



On the other hand, for an eigenmode with a complex 
eigenvalue, we rewrite the eigenmode as u = u' + iu" 
and v = v' + iv", where u',u",v' and v" are real, and 
we obtain 



tan 9 = 



K + v[ - vU - vL,) 2 + « + < - "-1 - v'-if 
X (u[ - v[ vU) 2 + « - v'l -v!L x + v'U? ' 

(67) 

Equation (|57|) does not depend on the overall phase of the 
eigenmode, i.e., it is invariant under (u,v) — > e tx (u,v) 
for an arbitrary real \- 

For the case of k = 0, 9 for the (u x , v x ) mode is analyt- 
ically obtained both for real and imaginary eigenvalues 
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Then the initial state is given as 



tan 2 6> 



1 - go 
1 + go 

1 - go 
1 + go 



fk - 9k 



Ik + 9k 



|^ +c dd n(l + go)(Q fc - Co) 



2M 



+ g-c dd n(l-g )(2Cfe + Co) 



(68) 



We also calculate numerically for k^O. In Fig. [31 we 
plot 9 for the lowest-energy modes which become unsta- 
ble in the red regions of Fig. [2] 

The Bogoliubov analysis shown in Figs. [2] and [3J sug- 
gests that for large q (> -Edd) the MDDI favors a helical 
spin structure, since the dynamical instability that exists 
at A = oo is suppressed as A decreases, as shown in Fig. [5] 
Actually, the ground-state spin structure at large q is a 
spin helix as discussed in Sec. IVI CI On the other hand, 
for small q (< E^d), fluctuations in both longitudinal and 
transverse magnetizations become dynamically unstable. 
The Landau instability at A = oo reflects the fact that 
E K has the minimum for k z ^= 0: 8 corresponding to the 
Landau instability is small [Figs. [3] (j) and (k)], imply- 
ing that the instability is caused by the fluctuations of 
transverse magnetizations. 

Although the Bogoliubov analysis indicates that the 
MDDI favors nonuniform spin structures.it is insufficient 
to account for the Berkeley experiment 0]: as discussed 
in Rcf. [6], the minimum wave length of the unstable 
modes (~ 30 /zm) is about three times larger than that 
of the observed magnetic pattern (A CX p ~ 10 /um). 



IV. SPIN DYNAMICS IN A TRAPPED SYSTEM 

To take into account the effects of the trapping poten- 
tial, initial conditions, noises, and nonlinearity, we per- 
form numerical simulations of the full 3D GPE with the 
dipole kernel in Eq. (|32p . We first consider a large ax- 
isymmctric pancake-shaped BEC to eliminate the effect 
of anisotropy of the trap, and then discuss the case of the 
Berkeley experiment [2| in the next section. 

We consider a BEC of N = 1.0 x 10 7 atoms 
in a harmonic trap with frequencies (yxiVyi v z) = 
(4.2,420,4.2) Hz. The corresponding peak density is 
n(0) = 2.0 x 10 14 cm" 3 and the Thomas-Fermi (TF) radii 
are (r x ,r y ,r z ) = (140,1.4,140) //m. The width along 
the y direction is smaller than the spin healing length 
£ sp = h/ y/2M\c\\n{0) = 2.8 /im and the dipole healing 
length £dd = V \/^Mcddn(0) = 9.5 //m. The magnetic 
field B = 165 mG is applied in the z direction. 

In the simulation, we first calculate a stationary state, 
^r( lm )(y) ! polarized in the m = — 1 state by using the 
imaginary-time propagation method, i.e., we have solved 
the stationary state of Eq. (|12p by replacing t with —it. 




N~\¥ ini) (r)\- 



(69) 



where n = kz, q K is defined in Eq. (|49|) whose denom- 
inator is replaced by 2n(0)(|cx| — CddQre), and N is a 
normalization constant; the parameters a, f3 and 7 are in- 
troduced to simulate fluctuations in the transverse mag- 
netization, the longitudinal magnetization, and the over- 
all phase, respectively, due to the quantum and thermal 
fluctuations as well as experimental noises. They are as- 
sumed to take on real random numbers independently 
on each grid and to obey the Gaussian distribution with 
variance cr Qj( 3. 7 . In the following calculation, we choose 
cr a = <7y = 0.1 and a/3 = 0.03. The initial noise depen- 



dence of the spin dynamics is discussed in Sec. IVI 



A. Effect of the nonuniform density 

We first consider the spin dynamics without the MDDI. 
Interestingly, a periodic pattern of the transverse magne- 
tization develops even in the case of Cdd = and k = 0. 
Figure S] shows the time evolution of the magnetic struc- 
ture for q/h = 2.9 Hz. Shown are (a) the direction of the 
transverse magnetization arg(/ + ), (b) longitudinal mag- 
netization f z /n(0), and (c) spin correlation function 



9(r±) 



= /d 2 rl/ + (rl+r ± )/_(rl) 
/ d 2 r' ± n(r' ± + r±)n(r'j_) 



(70) 



where n(r±) = J dyn(x,y, z) is the column density and 
f( r ±) = J dyf(x,y,z) is the column spin density. As 
shown in Fig. |4l the magnetic pattern develops from a 
uniform spin structure, and the periodic pattern appears 
in the spin correlation function at t = 1.8 s, which is 
destroyed in the further time evolution. This instability 
is due to the nonuniform density profile in the trapping 
potential. 

In the present case of a 87 Rb BEC, the density-density 
interaction dominates the system (cqh \c\\n, g), and 
therefore the density distribution n(r) is determined in- 
dependently of the spin configuration. Using Eq. (|36p . 
the ground-state order parameter in the local density ap- 
proximation (LDA) is given by 



*LDA(r) 



( 



2\ C1 \n(r) 



\ 



1 4- 1 

i+ 2| Cl |n(r) 



(71) 



for n(r) > g/(2|ci|), and 



LDA 



(r) = \/n(r) 



(72) 
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(a) arg(/+) 
t=0.72 sec 



(b) /*/n(0) 
*=0.72 sec 



t=1.8 sec 



*=3.2 sec 




t=1.8 sec 



t=3.2 sec 



-7T/3 




FIG. 4: (Color) Time evolution of (a) the transverse magne- 
tization, (b) the longitudinal magnetization, and (c) the spin 
correlation function (|70[) m the absence of the MDDI. The 
initial spin configuration is uniform and q/h — 2.9 Hz. The 
solid circles in (a) and (b) represent the TF radius (140 /im) 
at y — 0. The size of each panel is 310 /iin x 310 /im. The 
peak-to-peak distance of the correlation function at t — 1.8 s 
is about 70 /im. 



for n(r) < g/(2|ci|). The transverse magnetization per 
particle for \&lda distributes according to 



B. Effect of spin current 

The effect of the trapping potential becomes more 
prominent when the initial condition is a spin helix. Fig- 
ure [5] shows the spin dynamics in the absence of the 
MDDI for q/h = 1.8 Hz starting with a spin helix with 
A = 2tt/k = 60 /im, where Figs. [5] (a) and (b) are the 
snapshots of the transverse and longitudinal magnetiza- 
tions, respectively, and Fig. [5] (c) shows the time evolu- 
tion of the amplitude of the longitudinal magnetization 
M z = J dr\f z (r)\/N. Note here that the helical spin 
structure induces a spin current of the longitudinal mag- 
netization defined by 



•spin 



2Mi 



(75) 



The initial spin helix induces the spin current of j z pm (t = 
0) = —hnnz/(2M), and hence, the longitudinal magne- 
tization is accumulated at the top and bottom of the 
condensate as shown in Figs. [5] (b) and (c). Moreover, 
the spin current is reflected at the edge of the conden- 
sate, generating an interference pattern as shown in the 
snapshots at t = 0.76 s in Figs. [5] (a) and (b). In the 
time evolution, the spin helix unwinds [Fig. [5] (a) 1.5 s] 
and winds again in the opposite sense [Fig. [5] (a) 3.4 s]. 
In this dynamics, the direction of the spin current is also 
inverted. In the long-time scale of a few seconds, the 
spin texture oscillates between helix and anti-helix con- 
figurations, leading to the oscillations of M z as shown in 
Fig. [5] (c). The local longitudinal magnetization can be- 
come larger for the smaller quadratic Zeeman energy, and 
therefore, the period of the oscillations becomes longer. 



!/+(*•)! 

n(r) 



2\ Cl \n(r) 







for n (r) > (?3) 
for n (r) < 



while that of the initial state (|6"9"|) is constant (except for 
the contribution from the noise term): 



\U(r)\ 
n(r) 



= \ 1- 



2|ci|n(0) 



(74) 



Therefore, the initial sate (f6"9")l is not stationary at the low 
density region, and the instability grows from the periph- 
ery [Fig. 2] (a) 0.72 s]. In the course of time evolution, the 
amplitude of the local magnetization oscillates at the pe- 
riphery and then the fluctuations begin to penetrate into 
the central region at t ~ 1.8 s. These fluctuations induce 
a periodic pattern in the correlation function as shown 
in the snapshot at t = 1.8 s in Fig. [4] (c), which lasts 
for about 1 s and eventually dissolves [Fig. H] (c) 3.2 s]. 
The peak-to-peak distance in the correlation function at 
t = 1.8 s is about 70 fim. 



C. Effect of MDDI 

Next, we consider the effect of the MDDI. Figure [5] 
shows the result for k = and q/h = 0.46 Hz, where 
Figs. |5] (a), (b), and (c) are the snapshots of trans- 
verse magnetization arg(/ + ), longitudinal magnetization 
f z /n(0), and amplitude of the magnetization |/|/n(0), 
respectively, and Fig. [5] (d) shows the time evolution of 
the amplitude of the transverse magnetization per parti- 



cle M± = J dr^Jf% + fy/N and that of the longitudinal 

magnetization M z . The Bogoliubov spectrum shown in 
Figs. [2] (j) and[3](j) predicts that the spin wave along the 
x direction grows with the most unstable wave length of 
41 fan in the time scale of /i/|Imefc| ~ 0.74 s. The di- 
rection of the magnetic fluctuation for the most unstable 
mode is tan# = 0.96, i.e., both the transverse and lon- 
gitudinal fluctuations grow. The numerical result shown 
in Fig. [6] agrees well with the Bogoliubov analysis: the 
stripe pattern of the longitudinal magnetization grows 
over the time scale of ~ 1 s where the domain size is 
20 /im. The fluctuation in the transverse magnetization 
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(a) arg(/+) 
4=0 sec 



£=0.76 sec t=1.5 sec 



4=3.4 sec 




(a) arg(/+) 

4=0.76 sec 4=1.5 sec 4=3.4 sec 4=5.3 sec 



(b) 7z/n(0) 

4=0 sec 4=0.76 sec 4=1.5 sec 4=3.4 sec 



(c) M 




FIG. 5: (Color) Spin dynamics in the absence of the MDDI 
starting with a spin helix of A = 60 /im. Shown are the 
snapshots of (a) the transverse magnetization and (b) the 
longitudinal magnetization for q/h = 1.8 Hz, and (c) the time 
evolution of M z for q/h = 0.46 Hz (red), 1.8 Hz (green), 
and 2.9 Hz (blue). The size of each panel in (a) and (b) is 
310 /an x 310 /im. 



also grows. As the magnetic domain of the longitudinal 
magnetization develops, the magnetization |/| decreases 
at the domain walls [Fig. [6] (c) 1.5 s]. In the course of 
time evolution, these domain wall are destroyed by gener- 
ating pairs of polar-core vortices. The polar-core vortex 
is a topologically stable spin vortex whose core is un- 
magnetized, i.e., filled with the polar state. Each blue 
dot in the snapshots at t = 3.4 s and 5.3 s in Fig. [5] (c) 
indicates the core of a polar-core vortex. The obtained 
stripe structure is relatively stable with the domain size 
gradually becoming larger as shown in the snapshot at 
t = 5.3 s in Fig. [51(b). 

On the other hand, when q/h = 2.9 Hz and k = 0, the 
spin-wave mode of the transverse magnetization along 
the z direction becomes unstable [Figs. [5] (1) and[3](l)]. 
Figure[7]shows the numerical results for q/h = 2.9 Hz and 
K = 0, in agreement with the Bogoliubov analysis. Due 
to the large quadratic Zeeman energy, the longitudinal 
magnetization cannot grow for q/h = 2.9 Hz [Fig. 0(c)]. 

The initial spin helix changes the distribution of the 
unstable mode in the Bogoliubov spectrum (Fig. [2]), as 
well as induces the spin current of the longitudinal mag- 
netization. Figure [3] shows the spin dynamics starting 
with a spin helix with (a), (b) A = 120 /tm and with (c), 




(b) him 

4=0.76 sec 4=1.5 sec 4=3.4 sec 4=5.3 sec 




(c) |/|/fiC0) 

4=0.76 sec 4=1.5 sec 4=3.4 sec 4=5.3 sec 



1.0 



0.3 



(d) 




FIG. 6: (Color) Spin dynamics in the presence of the MDDI 
for q/h = 0.46 Hz starting with a uniform spin structure. 
The corresponding unstable mode is shown in Figs. [2] (j) and 
[3] (j). Shown are snapshots of (a) transverse magnetization, 
(b) longitudinal magnetization, and (c) amplitude of the mag- 
netization, and (d) time developments of M± (red) and M z 
(green). The size of each panel is 310 /im x 310 /im in (a) and 
(b), and 155 /im x 155 /im in (c). 



(d) A = 60 fim for q/h = 0.46 Hz. For A = 120 /im, the 
instability along the x direction grows. At the same time, 
the longitudinal magnetization is accumulated at the top 
and bottom of the condensate, leading to the magnetic 
pattern shown in the snapshots at t = 1.5 s in Figs. [8] (a) 
and (b). Then, helical structure is completely destroyed. 
On the other hand, for A = 60 /mi, the large initial spin 
current dominates the initial spin dynamics. Different 
from the case of Fig. [3J the pitch of the helix becomes 
smaller and smaller, and finally the helical structure is de- 
stroyed by generating pairs of polar-core vortices [Fig. [8] 
(c) 1.5 s]. In both cases of A = 120 /im and 60 /tm, the 
local longitudinal magnetization substantially increases 
[Fig. [5] (e)] , and the magnetic domains tend to elongate 
in the z direction. The polar-core vortices are located 
mainly at the domain wall of the longitudinal magneti- 
zation for both A = 60 /tm and 120 /im. We have also ob- 
served Mermin-Ho (MH) vortices for A = 120 /tm, where 
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FIG. 7: (Color) Spin dynamics in the presence of the MDDI 
for q/h — 2.9 Hz starting with a uniform spin structure. The 
corresponding unstable mode is shown in Figs. [2] (1) and [3] 
(1). Shown are snapshots of (a) transverse magnetization and 
(b) longitudinal magnetization, and (c) time evolution of M± 
(red) and M z (green). The size of each panel in (a) and (b) 
is 310 /jm x 310 /xm. 



the core of the MH vortex is magnetized and the direc- 
tion of the transverse magnetization changes ±2tt around 
it. The examples of the MH and polar-core vortices are 
indicated in Fig. [8j 



Figure [9] shows the spin dynamics for q/h = 2.9 Hz 
starting from a spin helix with a pitch of (a), (b) A = 
120 /xm and (c), (d) A = 60 /xm. Although the BEC for 
these parameters is stable in an infinite quasi-2D system 
[Fig. 0(f) and (i)], magnetic structures evolve due to the 
effect of the nonuniform density and the spin current. In 
the case for a large helical pitch of A = 120 /xm, the heli- 
cal structure of transverse magnetization is stable [Fig. [9] 
(a)]. In this case, a regular magnetic pattern of the longi- 
tudinal magnetization emerges as an interference pattern 
of the spin current [Fig. 0(b)]. On the other hand, in the 
case of A = 60 /xm, the large initial spin current induces 
the untwisting and rc-twisting of the helix as in the case 
of Fig. [5j although this dynamics is no longer periodic 
[Fig. 1(e)]. 



(a) A=120 /jm, arg(/+) 

£=0.76 sec £=1.5 sec t= 2.3 sec £=5.3 sec 




(b) A=120 ttm, f z /n(0) 

t=0.76 sec i=1.5 sec £=2.3 sec /" t=5.3 sec 





11.0 



-1.0 



(c) A=60 IMS, arg(/+) 

i=0.76 sec £=1.5 sec £=2.3 sec £=5.3 sec 




(d) A=60 Aim, f z /n(0) 

£=0.76 sec £=1.5 sec £=2.3 sec / £=5.3 sec 






(e) 



1.0 
0.8 
0.6 
0.4 
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M ± (A=120 /mi) 




M z (A=120 urn) 




Ml (A=60 ;im) 




M z (A=80 imi) - 



2 3 
time [s] 



FIG. 8: (Color) Spin dynamics in the presence of the MDDI 
for q/h = 0.46 Hz starting with a spin helix with a pitch of (a), 
(b) A = 120 £im and (c), (d) A = 60 /xm. Shown are snapshots 
of (a), (c) transverse magnetization and (b), (d) longitudinal 
magnetization, (e) Time evolution of M± and M z . The size 
of each panel in (a)-(d) is 310 /xm x 310 /xm. Examples for the 
MH vortex (MHV) and polar-core vortex (PCV) are enlarged 
at t — 2.3 s in (a), (b) and (c), (d), respectively, where the 
vortex core of MH vortex is magnetized and that of polar-core 
vortex is un-magnetized. In both cases of MH and polar-core 
vortices, the transverse magnetization changes 2ir around the 
vortex core. 
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(c) A=60 inn, arg(/ + ) 
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FIG. 9: (Color) Spin dynamics in the presence of the MDDI 
for q/h = 2.9 Hz starting with a spin helix with a pitch of (a), 
(b) A = 120 £im and (c), (d) A = 60 /jm. Shown are snapshots 
of (a), (c) transverse magnetization and (b), (d) longitudinal 
magnetization, (e) Time evolution of M± and M z . The size 
of each panel in (a)-(d) is 310 /im x 310 /im. The sense of 
twisting of the spin helix dose not change in (a), while it is 
reversed in (c) at t = 2.3 s. 



COMPARISON WITH THE BERKELEY 
EXPERIMENT 



A. Numerical results 

Now we consider the system of the Berkeley experi- 
ment Q. The difference from the previous subsection 
is that (i) the trapping potential is elliptical and elon- 
gated along the direction of the magnetic field, (ii) the 
initial state is fully magnetized [|/(r)| = n(r)} while 



tcr (|69|) . and (iii) the spin helix is evolved by applying a 
magnetic field gradient during a period of 5-8 ms. 

In accordance with the Berkeley experiment, we con- 
sider a BEC of TV = 2.3 x 10 6 atoms in a harmonic 
trap with frequencies (u x , u y , u z ) = (39, 440, 4.2) Hz. The 
peak density is n(0) = 2.8 x 10 14 cm~ 3 and the TF radii 
are (r x ,r y ,r z ) = (18, 1.6, 169) fim. The spin and dipolc 
healing lengths arc £ sp = 2.3 /im and £dd = 8.0 [im, re- 
spectively. 

For the initial state, we prepare a spin-polarized state 
in the x direction as 




= JV|#0 ni >(r)|e* 



r ia 1-P 
c 2 




5 - 4, 



(76) 



|/(r)| 



< n(r) for the order parame- 



where \I/( ini )(r), Af, a, (3, and 7 are the same as those ap- 
pearing in Eq. (|6"9"|) . We introduce S to simulate fluctu- 
ations in the amplitude of the magnetization. Here we 
assume that S takes on a complex random number inde- 
pendently on each grid and obeys the Gaussian distribu- 
tion with variance as- In the following numerical calcu- 
lation, we choose <j a = cr 7 = 0.1 and op = 75 = 0.03. 
We prepare a helical structure with pitch A in the real- 
time evolution by applying a field gradient of dB/dz = 
h/ (\gp\fiB^Tp) during a period of t p = 5 ms. 

We first consider the spin dynamics for n = and 
Cdd = 0. Since the TF radius in the x direction is 
small, the effect of the nonuniform density becomes more 
prominent than the previous case. In addition, the 
deviation from the LDA stationary state discussed in 
Eqs. ([7T |) - (f74| becomes larger since we start with a fully- 
magnetized state. Hence, the time and length scales of a 
checkerboard pattern become smaller; the periodic pat- 
tern emerges spontaneously at t ~ 0.7 s and the peak-to- 
peak distance of the correlation function is about 25 /im. 
The snapshots of the spin configuration at t = 0.82 s are 
shown in Fig. [TU] (a). The periodic pattern eventually 
dissolves in a few seconds. 

When we start from a spin helix with Cdd = 0, the lon- 
gitudinal magnetization is accumulated at the top and 
bottom of the condensate as in the case of a pancake- 
shaped trap. However, different from the pancake-shaped 
trap, the checkerboard pattern of the transverse magne- 
tization with the domain size 15 ~ 20 fim appears at the 
top and bottom edges of the condensate for t > 0.7 s due 
to the strong effects of the nonuniform density and the 
interference of the spin current. The oscillations between 
the helix and anti-helix structures are observed as in the 
case of the pancake-shaped trap (Fig. [5]), and the checker- 
board pattern disappears during these oscillations. 

Figure QT] shows the spin dynamics in the presence of 
the MDDI for k = 0. For q/h = 0.46 Hz, fluctuations 
of the longitudinal magnetization grow in the x direction 
[Fig. [TT] (a) 0.82 s], in agreement with the Bogoliubov 
analysis [Figs. [5] (j) and|3](j)]. However, in the further 
time evolution, the domains are aligned along the z direc- 
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(a) no helix, t 
arg(7+) 7z 



=0.82 sec (b) helix. 

sec 



A=60 fj,m 
0.77 1.3 2.0 




-7T _ I 7T 

arg(/+) 

FIG. 10: (Color) (a) Snapshots of spin configuration at 
t = 0.82 s in the absence of the MDDI starting with a uniform 
spin configuration for q/h = 2.9 Hz, where each panel shows 
the distribution of transverse magnetization arg(/+) (left), 
longitudinal magnetization f z /n to t(0) (middle), and spin cor- 
relation function g (right). The peak-to-peak distance of the 
correlation function is about 25 fim. (b) Time development of 
transverse magnetization in the absence of the MDDI start- 
ing with a spin helix of pitch A = 60 fim for q/h — 1.8 Hz. 
The domain size of the checkerboard pattern in (b) is about 
20 /im. The solid curves in the distribution of longitudinal 
and transverse magnetizations represent the TF boundaries 
at y = 0. The size of each panel is 45 /xm x 359 /xm. Note 
that the color scales of arg(/+) in (a) and (b) are different. 



tion [Fig.[TT](a) 3.1 s], because the long domain wall par- 
allel to the z direction is energetically unfavorable. This 
domain structure seems to be a micro-canonical equilib- 
rium state for small q, and lasts for a long time [Fig. QT] 
(a) 5.1 s]. In the equilibrium state, the transverse mag- 
netization forms a helix pattern along the z direction 
[Figs. QT] (b)]. At large q, the checkerboard pattern first 
develops in the short time scale of t < 1 s as in the case 
of Cdd = 0. Then, the effect of the MDDI on the spin 
dynamics becomes more prominent in the longer time 
scale. For q/h = 2.9 Hz, the fluctuations of magnetiza- 
tions grow mainly in the transverse direction and form 
a helix along the z direction after a few seconds [Fig. [TT] 
(c)]- 

Finally, we show the result for the case of the Berkeley 
experiment Q, i.e., the spin dynamics in the presence of 
the MDDI starting with a spin helix. Figure [T^] shows the 
result for A = 60 /im and q/h = 1.8 Hz. In the short time 
scale of t < 1 s, the longitudinal magnetization is accu- 
mulated at the top and bottom of the condensate due to 
the spin current. At the same time, the helical structure 
is modulated at the top and bottom of the condensate 
due to the effects of the nonuniform density and intcrfcr- 



(a) g=0.46 Hz, f z 
0.82 sec 2.0 2.7 



3.1 



5.1 



(c) 
9=2.9 Hz 
t=3.3 sec 
fz arg(j>) 



fz/n(0) 



FIG. 11: (Color) Magnetic structures developed from a uni- 
form spin configuration in the presence of the MDDI at (a), 
(b) q/h = 0.46 Hz and (c) 2.9 Hz, where (a) and the left panel 
in (c) show the distributions of the longitudinal magnetiza- 
tion, while (b) and the right panel in (c) show the direction 
of the transverse magnetization. The longitudinal magnetiza- 
tion is scaled with /f ax = 1.0 in (a), and ]f ax = 0.6 in the 
left panel of (c). The transverse magnetization forms a helical 
structure in both cases for q/h = 0.46 and 2.9 Hz. The size 
of each panel is 45 /jm x 359 /tin. 



ence of the spin current. The pitch of the helix becomes 
smaller and smaller, and the helix is eventually destroyed 
due to the MDDI. In this dynamics, both polar-core vor- 
tices and MH vortices are generated spontaneously as 
indicated in Fig. [T2j 



B. Comparison with the experiment 

Here we summarize the agreements and disagreements 
between the Berkeley experiment 0] and our numerical 
simulation. 

In agreement with the experiment, the uniform spin 
structure is stable over a few hundreds of milliseconds. 
However, unlike the experimental results, the spin helix 
is also stable within the time scale of a few hundreds of 
milliseconds with and without the MDDI. In the longer 
time scale of a few seconds, magnetic patterns develop 
from both the uniform spin structure and the spin helix 
due to the trapping potential and the MDDI. Even in 
the absence of the MDDI, the nonuniform density dis- 
tribution and interference of the spin current induce the 
checkerboard pattern (Fig. [TU)) . However, the domain 
size observed in this dynamics is at least three times 
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FIG. 12: (Color) (a) Time development of the direction of the 
transverse magnetization with the MDDI starting with a spin 
helix of pitch A = 60 fim for q/h = 1.8 Hz. (b) Snapshot of 
the longitudinal magnetization at t — 4.0 s. The size of each 
panel is 45 /im x 359 fim. The solid circles and dotted circles 
in (a) 4.0 s and (b) indicate the locations of the polar-core 
vortex (PCV) and MH vortex (MHV) , respectively. 



larger than that observed in the experiment, and the 
checkerboard pattern eventually dissolves in our simula- 
tions. The MDDI does not stabilize this periodic pattern. 
The pattern induced by the MDDI is sensitive to the 
quadratic Zeeman energy. When the MDDI dominates 
the quadratic Zeeman energy, the staggered domain of 
the longitudinal magnetization appears [Figs.[TT](a)]. On 
the other hand, when the quadratic Zeeman energy dom- 
inates the MDDI, magnetization is almost transverse and 
forms a spin helix [Figs.Qj](c)]. The length scale of these 
structures is of the order of some tens of micrometers and 
much larger than that observed in the experiment. 

Apart from the time and length scales, there are several 
discrepancies in the property of the magnetic structure 
between the experiment and our calculation. First, in 
our calculation, the instability is always accompanied by 
the emergence of the local longitudinal magnetization, 
whereas in the experiment the longitudinal magnetiza- 
tion is much smaller than the transverse one. In partic- 
ular, when we start with a spin helix, the longitudinal 
magnetization grows rapidly and becomes comparable to 
the transverse magnetization. Second, the spin structure 
induced by the MDDI is sensitive to the external mag- 
netic field, whereas the spin dynamics observed in the 
Berkeley experiment is insensitive to the quadratic Zee- 
man energy for 0.8 < q/h < 4 Hz. Third, the growth 
rate of the checkerboard pattern (due to the nonuniform 
density and the spin current in our case) is insensitive 
to the pitch of the initial helix, whereas the growth rate 



increases as the pitch becomes smaller in the experiment. 
On the other hand, in agreement with the experiment we 
have observed the spontaneous generation of pairs of the 
polar-core vortices. 



VI. DISCUSSION 

We here consider the possible reasons for the discrep- 
ancies discussed in the previous section. 



A. Initial noise dependence 

To simulate fluctuations and noises, we have intro- 
duced the initial noises of a, /3, 7 and 6 in Eqs. (1691) and 
(|76l) . Here we calculate the spin dynamics for various 
variances of initial noises. 

When we start with a uniform spin structure, the fluc- 
tuations in the transverse (a) and longitudinal (/3) mag- 
netizations contribute to the growth of the instability, 
while the fluctuations in the overall phase (7) and am- 
plitude of magnetization (5) hardly affect the formation 
of spin textures. In Figs. Q2] (a) and (b), we show the 
time evolution of the magnetization starting from a uni- 
form spin structure with various initial noises. The fluc- 
tuation does not grow for a a = ap = during a few 
seconds. This result is consistent with the fact that the 
unstable mode in an infinite quasi-2D system is propor- 
tional to (1, 0, — 1) T and decouples from the fluctuations 
in the overall phase (phonon) and the amplitude of mag- 
netization. Due to the effect of the trapping potential, 
the fluctuation eventually grows after a few seconds for 
a a = 073 = 0. 

On the other hand, when we start from a spin helix, 
the spin dynamics is independent of the detail of the 
initial noise, since the spin current dominates the initial 
dynamics [Figs. Q2] (c) and (d)]. 

Although the magnetic pattern slightly depends on the 
initial noise, the dependence is too little to account for 
the discrepancy between the experiment and the numer- 
ical result discussed in the previous section. Figure [T4l 
shows an example of the initial noise dependence on the 
spin textures, where the magnetic patterns at t = 5.3 s 
with the initial noise of (a) a a — 0.03 and erg = er 7 = 
as = and (b) ap — 0.03 and a a = <r 7 = as = are 
shown. In both cases, M z saturates to ~ 0.7 and similar 
magnetic structures develop. However, when the initial 
noise is transverse (longitudinal), the length scale of the 
spatial structure of the transverse (longitudinal) magne- 
tization is smaller than the case when the initial noise is 
longitudinal (transverse), which means the kinetic energy 
is stored in the fluctuation of the transverse (longitudi- 
nal) magnetization. 
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FIG. 13: (Color online) Time evolution of spin fluctua- 
tions of (a), (c), (d) M z = J dr\f z (r)\/N and (b) M y = 
| J drf y (r)\/N , where a a — 0.1 means <r a = 0.1 and other 
components vanish, etc. In (c) and (d), a = and 0.1 mean 
o~a = cp = o 1 — as = and 0.1, respectively. 



FIG. 14: (Color) Magnetic structure at t = 5.3 s starting 
from a uniform spin structure for q/h — 0.46 Hz with the 
initial noise of (a) o a — 0.03 and ap = a 7 = = and (b) 
aa = 0.03 and <j a — <r 7 = ag — 0. The length scale of the 
fluctuation in the transverse magnetization is smaller (larger) 
than that in the longitudinal magnetization in (a) [(b)]. The 
size of each panel is 45 /im x 359 fim. 



B. Uncertainty of the parameters 

We examine the effect of the uncertainty of the coef- 
ficient c\ . The strength of the spin exchange interaction 
has been measured in molecular spectroscopy [29j j and in 
spin-mixing dynamics (3ll . |32| . According to these mea- 
surements, spin-exchange interaction energy \c\\n lies be- 
tween 7 Hz and 13 Hz for n = 2.8 x 10 14 cm~ 3 , where 
the error bar in the difference of the scattering lengths 
is — l.OaB < a-2 — ao < — 1.8aB- A decrease in |ci| en- 
hances the quadratic Zeeman effect. Therefore, the effect 
of a nonuniform density distribution is enhanced and the 
longitudinal magnetization induced by spin helix and the 
effect of the MDDI are suppressed. However, the prop- 
erty of the magnetic pattern is qualitatively unchanged 
within the above range of the uncertainty and it cannot 
resolve the discrepancy between the Berkeley experiment 
and our calculation. 

The uncertainty of the uniform magnetic field does 
not affect the spin dynamics as long as the magnetic 
field is much stronger than the dipole field, since we 
use the time-averaged dipole kernel (|33[) in the numer- 
ical simulations. On the other hand, the residual field 
gradient may change the dynamics. We calculate the 
spin dynamics under the residual field gradient dB/dz = 
10 [iG/(2ztf) [HI]. However, the property of the mag- 
netic pattern is almost unchanged in the short time scale 
of a few hundreds milliseconds, although the residual field 
gradient winds helix more and more in the long time 
scale. 



C. Stable spin structure with the MDDI 

Next we consider the effect of the energy dissipation on 
the dynamics by replacing t with (1 — iT)t in Eq. (|12l) . 
To keep the total longitudinal magnetization and number 
of atoms constant, we introduce the term pm^> m in the 
right-hand side of Eq. (| 12[) with p being the Lagrange 
multiplier, and change p and the chemical potential fi in 
each step. The energy dissipation leads to enlargement 
of the spatial structure rather than stabilization of the 
patterns, regardless of the value of T over the range of 
0.001 < T < 0.1, and eventually the system reaches a 
stationary state with a spatial structure of the order of 
100 /im. 

We choose the same trap frequencies (lj x , u) y , uj z ) = 
2tt x (39,440,4.2) Hz and the number of atoms N = 
2.3 x 10 6 as used in Ref. 0], and investigate the stable 
spin configuration for various q. The obtained results 
are shown in Fig. 1151 and the q dependence of the ampli- 
tude of the transverse, longitudinal, and total magnetiza- 
tions are plotted in Fig. [TBI We have found two types of 
stable spin textures, and the first-order phase transition 
between these structures occurs at q = q c = h x 2.3 Hz. 
The critical quadratic Zeeman energy q c is close to MDDI 
energy E^d — Cdd^So; for the present system with 
n = 2.8 x 10 14 cm~ 3 , E dd ~ h x 1.9 Hz. 

The spin configuration is determined by the interplay 
between the MDDI and the quadratic Zeeman effect. 
When the MDDI dominates the quadratic Zeeman en- 
ergy, the longitudinal magnetization is favored since the 
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spin-dcpcndcnt factor 8 vv i — i5 vz 5 v i z in Eq. (|32|1 con- 
tributes maximally for v = v' = z. From the orbital part 
in Eq. (|32|) . we find that the magnetic domain of the lon- 
gitudinal magnetization elongates in the direction of the 
external magnetic field. The condensate is almost fully 
magnetized and the magnetization at the domain wall is 
perpendicular to the magnetic field. The direction of the 
domain wall is determined by the kinetic energy of the 
domain wall, and depends on the aspect ratio of the trap: 
the domain wall is perpendicular to the magnetic field for 
the trap used in the Berkeley experiment whereas it is 
parallel to the magnetic field in a pancake-shaped trap. 

As the quadratic Zeeman energy becomes larger, the 
width of the domain wall becomes larger as shown in 
Fig. [TS] (a), where the magnetization at the domain wall 
and at the top and bottom of the condensate is perpen- 
dicular to the magnetic field. For q > q c , the quadratic 
Zeeman energy dominates the MDDI, and magnetiza- 
tion occurs in the plane perpendicular to the magnetic 
field. As regards transverse magnetization, the MDDI 
favors the antiferromagnetic ordering in the direction of 
the magnetic field, resulting in the helical spin structure 
as shown for the case of q/h = 2.6 Hz in Fig. [15] (b). 

The ground-state spin structure has also been inves- 
tigated in Ref. Q in a quasi- 2D system by using the 
Metropolis Monte Carlo method. Our results in Figs. [15] 
and [TH] are consistent with those in Ref. Q • 
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FIG. 15: (Color) Stationary configuration of (a) longitudinal 
magnetization / 2 /ntot(0) and (b) transverse magnetizations 
arg(/+) for the trap geometry and the number of atoms used 
in Ref. 0. The condensate is almost fully magnetized in all 
cases (See Fig. ll6|l , The size of each panel is 45 /im x 359 (jm. 

We have also searched a metastable state with a pe- 
riodic pattern, including the vortex lattice state of both 
MH and polar-core vortices suggested in Ref. Q. We 
add a periodic modulation to the initial order parame- 
ter and investigate its dynamics with the imaginary-time 
propagation. However, all configurations that we have 
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FIG. 16: (Color online) q dependence of the amplitude of the 
transverse magnetization M±_ = f dr^J f% + fy/N, longitudi- 
nal magnetization M z — J dr\f z \/N, and total magnetization 
M = f dr-J fi + fy + fi/N, of the stationary state for the 
trap geometry and the number of atoms used in Ref. 0]. 



prepared were unstable and eventually goes to the struc- 
tures shown in Fig. 1151 



D. Spin vortex lattice 

It is pointed in Ref. Q that the vortex lattice of the MH 
vortices is long-lived in the imaginary-time propagation; 
unfortunately we have not been able to reproduce such 
tendency. To evaluate the lifetime of the vortex lattice, 
we here investigate the real-time dynamics starting with 
the vortex-lattice state. As mentioned before, there are 
two types of vortices in this system: the polar-core vortex 
and the MH vortex. The order parameter around a polar- 
core vortex is given by 



o±i<P 



PCV 



fv(r) 



' fv(r) 



(77) 



where we assume that the order parameter is axisymmet- 
ric around the vortex; (p is an azimuthal angle around 
the vortex, and /„ (r) is a monotonically increasing func- 
tion which satisfies /„(0) = and lim r _>oo f v (r) = 1/2. 
The spin current defined in Eq. (|75|) is anti-clockwise 
for 1 3>p cv and clockwise for \I>£ CV 5 so that they form 
a vortex-antivortex pair. The order parameter around a 
MH vortex is given by 



o ll P 



/ e ±i<p 1+cos ff±(r) 
sin /3± (r) 



MHV 



V2 

l-cosfi ± 0) 

„±itp 1+cos ff T (r) ' 
2 

sin/3 T (r) 
%/2 

„Ti<g l-cos^ T (r) 



(78) 



(79) 
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where /3+ (r) and /3_ (r) are monotonically decreasing and 
increasing function, respectively, which satisfy /3+(0) = 
7T, /3-(0) = and linv^oo (3± (r) = ir/2. The mass cur- 
rent 



2Mi 



; V* m - (V^)*J 



(80) 



is anti-clockwise for \I>mhv an( ^ c i 0C ]j W i se f or <j>mhv^ 

MHV and 



while the spin current is anti-clockwise for ty± , 
clockwise for $™ v . Figure [T7] shows the spin struc- 
ture around each vortex in Eqs. (|77jl (|79|) and the possi- 
ble configuration for the vortex lattice with periodic spin 
structure. In order for the spin configuration to be peri- 
odic, the vortices with clockwise and anti-clockwise spin 
current have to align alternately as shown in Fig. 1171 (c). 
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FIG. 17: (a) Spin configuration around vortices in Eqs. (|77[) - 
(|79|l . where the arrows show the direction of the transverse 
magnetization, P indicates the polar core, and Q and 
mean that the magnetization at the vortex core points +z 
and —z direction, respectively, (b) Two configurations are 
topologically equivalent if the core structure is the same, (c) 
Unit of a vortex lattice that exhibits a periodic structure of 
the transverse magnetization, (d)-(e) Possible sets of vortices 
to form a vortex lattice. 



we use the ansatz f v (r) = ^tanh(r/£ sp ), f3+(r) = n — 
■|tanh(r/^ sp ), and /3_(r) = -|tanh(r/^ sp ). Among the 
configurations in Figs. [17] (d)-(g), Fig. [17] (f) is the most 
unstable and dissolves within 200 ms, since this config- 
uration has a nonzero mass circulation. On the other 
hand, the periodic pattern survives longest for the con- 
figuration shown in Fig. [17] (e). Figure [18] shows the 
spin dynamics starting from the vortex lattice in Fig. [T7] 
(c) with lattice constant d v = 10 /im. In this case, the 
periodic pattern in the correlation function dissolves at 
around 0.5 s, both in the presence and absence of the 
MDDI. 

Our results show that if the MH vortex-lattice emerges 
for some reason, it can survive for 0.5 s, which is longer 
than the time scale for the emergence of the crystalline 
pattern in the experiment @]. However, since spin dy- 
namics in this system is slow, for instance, the helical spin 
configuration lasts more than 0.5 s as shown in Fig. [T2l 
we cannot conclude that the spin vortex lattice is long- 
lived compared with other configurations. Moreover, no 
mechanism has been presented for the vortex lattice to 
appear in a short time scale of a hundred milliseconds. 
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FIG. 18: (Color) Spin dynamics starting with spin-vortex lat- 
tice shown in Fig. [17] (e) at q/h = 1.8 Hz for the trap geometry 
and the number of atoms used in Ref. [§]. (a)-(c) Transverse 
magnetization arg(/_|_) (left), longitudinal magnetization f z 
(middle), and spin correlation function g (right) calculated in 
the presence of the MDDI. (d) Spin correlation function g cal- 
culated in the absence of the MDDI. Note that the periodic 
pattern in (d) is quite similar to that in the right panel of (b). 
The size of each panel is 45 /im x 359 /im. 



E. Thermalization due to the dipole relaxation 



We prepare the spin-vortex-lattice state shown in 
Figs. [T7] (d)-(g) with lattice constant d v , and investi- 
gate the spin dynamics in the real-time evolution. Here 



Finally, we discuss the validity of the time-averaged 
dipole kernel (|33|). It is pointed in Ref. [H that the 
terms in the MDDI which arc canceled by taking the 
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time average induce the instability of the Larmor preces- 
sion. The wave number of the unstable mode is about 
y/TMiJZ/h. This instability corresponds to the dipole 
relaxation, or the Einstein-de Haas effect fl9l - l2lj ]: when 
two atoms in the highest Zeeman sublevel (m = — 1) col- 
lide with each other, one or both of them can change 
their spin state to m = via the MDDI; due to the en- 
ergy conservation, the total kinetic energy of these atoms 
nearly equals to the linear Zeeman energy (Iiwl or 2Hujl 
depending on the number of spin- flipped atoms). The 
experimental result that the total longitudinal magne- 
tization is conserved for a period longer than the time 
scale of the MDDI [TtJ indicates that the effect of the 
dipole relaxation is quite small. However, since the lin- 
ear Zeeman energy hw^/k^ ~ 5.5 fiK is much higher 
than the temperature of the condensate (~ 100 nK) in 
the present system, a small number of spin-flipped atoms 
might thermalize the condensate. The treatment of the 
thermal atoms is beyond the Gross-Pitaevskii formalism, 
and remains a challenge for the future work. 

VII. CONCLUSIONS 

We have investigated the pattern formation dynam- 
ics in a spin-1 spinor dipolar Bose-Einstein condensate 
(BEC) observed by the Berkeley group 0, by taking 
into account the effects of spinor and dipolar interac- 
tions, linear and quadratic Zeeman energies, anisotropic 
trap geometries, various initial conditions, and noises. 

We have first performed the Bogoliubov analysis in a 
uniform quasi-two-dimensional system, and found that 
the spin helix does not enhance the instability but sta- 
bilize the system in the presence of the magnetic dipole- 
dipole interaction (MDDI) (Figs. [2] and ©. The mini- 
mum wavelength of the unstable mode is at least three 
times larger than the wavelength of the spin modulation 
observed in the experiment [2|. We have investigated 
the spin dynamics by simulating the three-dimensional 
Gross-Pitaevskii equation. There are three mechanisms 
that generate spatial spin structures: (i) a nonuniform 
density profile in an optical trap induces checkerboard 
pattern, even when we start from a uniform spin struc- 
ture in the absence of the MDDI [Fig. [10] (a)]; (ii) the 
initial spin helix induces a spin current, which is reflected 
at the edge of the condensate and generates a checker- 
boardpattern in the case for the trap geometry used in 
Ref. [2[ [Fig. [10] (b)], and (iii) in agreement with the Bo- 
goliubov analysis, the MDDI contributes to the pattern 
formation. However, in all cases, the domain size of the 
obtained magnetic pattern (~ 15 /xm at minimum) is 
more than three times larger than that observed in the 



experiment (A eX p/2 ~ 5 /mi). It takes more than 500 ms 
for the pattern to develop in our calculation, whereas it 
develops within 200 ms in the experiment. 

The other significant differences from the experiment 
are that (i) when the magnetic pattern develops, it is al- 
ways accompanied by the growth of the local longitudi- 
nal magnetization, and (ii) the MDDI-induced dynamics 
strongly depends on the strength of the quadratic Zee- 
man energy q. 

The detail of the initial noise does not qualitatively 
change the spin dynamics, as long as it includes the 
fluctuations in transverse or longitudinal magnetizations. 
We have also investigated the stationary spin structure in 
this system and obtained different spin textures from that 
observed in the Berkeley experiment. The stable texture 
undergoes the phase transition from a staggered domain 
of the longitudinal magnetization at small quadratic Zee- 
man energies to a spin helix of the transverse magne- 
tization at large quadratic Zeeman energies due to the 
interplay between the quadratic Zeeman energy and the 
MDDI energy. We have also considered the stability of 
the vortex-lattice state; although the checkerboard pat- 
tern can survive for about 500 ms, it is not long-lived 
compared with other spin textures. 

From the above discrepancies, we conclude that the 
mean-field and Bogoliubov theories at zero temperature 
cannot account for the Berkeley experiment 0- The ef- 
fects absent in our calculation are many-body correla- 
tions and thcrmalization. In particular, there might be a 
non-trivial effect of the thermal atoms via the dipole re- 
laxation as discussed in Sec. IVIEI Since the system size 
in the direction of the strong confinement is smaller than 
the spin healing length, the geometry of the condensate is 
two-dimensional with respect to the spin degrees of free- 
dom. Since the effects of fluctuations become prominent 
in low dimensional systems, thermal fluctuations might 
contribute significantly to the magnetic pattern. If ex- 
perimental external noises are not the origin of the mag- 
netic pattern, it would be a new quantum phenomenon 
beyond the mean-field theory. The clarification of these 
issues remains a challenge for a future work. 
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